%% Practical: Sensorspace Time-Frequency
%
% In this practical we will run through an analysis step by step, several
% of which require manual intervention. This will go through the following
% steps:
%
%   1) Prepare raw data for OAT analysis
%   2) Bandpass filter data and split into epochs
%   3) Compute a first level GLM analysis with OAT
%   4) Visualise results with FieldTrip
%   5) Compute a topoplot averaged within a time-frequency window
%
% We will work with a single subject's data from an
% emotional faces task (data courtesy of Susie Murphy) and perform an
% time-frequency analysis in sensor space. 
% This dataset can be downloaded from:  
%  
% www.fmrib.ox.ac.uk/~woolrich/faces_subject1_data.tar.gz
% 
% Note that this contains the spm file:
% spm8_meg1.mat
% that is an SPM MEEG object that has continuous data that has already been
% SSS Maxfiltered and downsampled to 250 Hz. 
% 
% and
% espm8_meg1.mat
%
% which is an SPM MEEG object that has the same data epoched into the
% different task conditions.
%

%% SETUP THE MATLAB PATHS
%
% Sets the Matlab paths to include OSL. Change these paths so that they 
% correspond to the setup on your computer. You will also need to ensure 
% that fieldtrip and spm are not in your matlab path (as they are included
% within OSL).

global OSLDIR;
    
% set this to where you have downloaded OSL:
osldir = '/Users/andrew/Software/Matlab/osl/osl_full/osl1.5.0_beta';  

addpath(osldir);
osl_startup(osldir);

%% INITIALISE GLOBAL SETTINGS FOR THIS ANALYSIS
%
% This cell sets the directory that OAT will work in. Change the workingdir variable to correspond to the correct directory on your computer before running the cell.

workingdir = '/Users/andrew/Projects/OSL_test/preproc_manual/faces_subject1_data';

cmd = ['mkdir ' workingdir]; if ~exist(workingdir, 'dir'), unix(cmd); end% make dir to put the results in

clear spm_files spm_files_epoched;

%% SET UP THE LIST OF SUBJECTS AND THEIR STRUCTURAL SCANS FOR THE ANALYSIS
%
% Specify a list of the fif files, structual files (not applicable for this
% practical) and SPM files (which will be created). It is important to make
% sure that the order of these lists is consistent across sessions. 
% Note that here we only have 1 subject, but more generally there would be
% more than one, e.g.:

% fif_files{1}=[testdir '/fifs/sub1_face_sss.fif']; 
% fif_files{2}=[testdir '/fifs/sub2_face_sss.fif']; 
% etc...

% spm_files{1} = [workingdir '/sub1_face_sss.mat'];
% spm_files{2} = [workingdir '/sub2_face_sss.mat'];
% etc...
% set up a list of SPM MEEG object file names (we only have one here)

spm_files{1}=[workingdir '/spm8_meg1.mat'];
spm_files_epoched{1}=[workingdir '/espm8_meg1.mat'];

%% SETUP SENSOR SPACE OAT SOURCE RECON
% 
% This stage sets up the source reconstruction stage of an OAT analysis.
% The source_recon stage is always run even for a sensorspace analysis,
% though in these cases it simply prepares the data for subsequent
% analysis.
%
% In this example we define our input files (D_continuous and D_epoched) 
% and conditions before setting a time frequency window from -400ms before
% stimulus onset to +500ms and 1 to 40Hz. The source recon method is set to
% 'none' as we are performing a sensorspace analysis

oat=[];
oat.source_recon.D_continuous=spm_files;
oat.source_recon.conditions={'Motorbike','Neutral face','Happy face','Fearful face'};
oat.source_recon.D_epoched=spm_files_epoched; % this is passed in so that the bad trials and bad channels can be read out
oat.source_recon.freq_range=[1 40]; % frequency range in Hz
oat.source_recon.time_range=[-0.4 0.5];
oat.source_recon.method='none';
oat.source_recon.dirname=[oat.source_recon.D_epoched{1} '_hilbert_multiband'];


%% SETUP SENSOR SPACE OAT FIRST LEVEL TF PARAMETERS
%
% Next we set up a single subject trial-wise GLM on our prepared data.
% Firstly the time-frequency parameters are defined, these must be within
% the bounds of the time-frequency window set in the source recon stage.
%

oat.first_level.tf_method='morlet'; % can be morlet or hilbert
oat.first_level.tf_freq_range=[5 40]; % frequency range in Hz
oat.first_level.time_range=[-0.2 0.3]; % need to make this time range smallet than oat.source_recon.time_range to remove edge effects
oat.first_level.tf_num_freqs=14; % we are keeping this unusally low in the practical for the sake of speed
%oat.first_level.tf_hilbert_freq_res=8;

% NOTE that you can also set HILBERT freq ranges explicitly, e.g.:
% oat.first_level.tf_hilbert_freq_ranges=[[4 8];[8 12];[12 16];[16 20];[20 24];[24 30]]; % frequency range in Hz

oat.first_level.post_tf_downsample_factor=4; % does downsampling after the TF decomposition

oat.first_level.bc=[1 1 0]; % specifies whether or not baseline correction is done for the different contrasts

%% SETUP SENSOR SPACE OAT FIRST LEVEL GLM
%
% This cell defines the GLM parameters for the first level analysis.
% Critically this includes the design matrix (in Xsummary) and contrast
% matrix
%
% Xsummary is a parsimonious description of the design matrix.
% It contains values Xsummary{reg,cond}, where reg is a regressor no. and cond
% is a condition no. This will be used (by expanding the conditions over
% trials) to croat_settingse the (num_regressors x num_trials) design matrix:
%
% Each contrast is a vector containing a weight per condition defining how
% the condition parameter estimates are to be compared. Each vector will
% produce a different t-map across the sensors. Contrasts 1 and 2 describe
% positive correlations between the each sensors activity and the presence
% of a motorbike or face stimulus respectively. Contrast 3 tests whether
% each sensors activity is larger for faces than motorbikes.

Xsummary={};
Xsummary{1}=[1 0 0 0];Xsummary{2}=[0 1 0 0];Xsummary{3}=[0 0 1 0];Xsummary{4}=[0 0 0 1];
oat.first_level.design_matrix_summary=Xsummary;

% contrasts to be calculated:
oat.first_level.contrast={};
oat.first_level.contrast{1}=[3 0 0 0]'; % motorbikes
oat.first_level.contrast{2}=[0 1 1 1]'; % faces
oat.first_level.contrast{3}=[-3 1 1 1]'; % faces-motorbikes
oat.first_level.contrast_name{1}='motorbikes';
oat.first_level.contrast_name{2}='faces';
oat.first_level.contrast_name{3}='faces-motorbikes';
oat.first_level.report.first_level_cons_to_do=[3 1 2];
oat.first_level.post_tf_downsample_factor=4;


%% RUN OAT
%
% Next we sanity-check our OAT options before running the source_recon and
% first_level stages of the analysis.

oat = osl_check_oat(oat);


oat.to_do=[1 1 0 0];
oat = osl_run_oat(oat);


%% LOAD IN RESULTS
%
% load first-level GLM result

stats=osl_load_oat_results(oat,oat.first_level.results_fnames{1});

% display freq bins used:
disp('Freq bins in Hz:');
disp(stats.frequencies);


%% VISUALISE USING FIELDTRIP
%
% Next we will use an osl wrapper around a Fieldtrip function to the
% results from contrast 3 (Faces>Motorbikes). This requires us to define
% several parameters in the S2 structure. Critically, the oat analysis and
% the number contrast within it.
%
% Note that this produces an interactive figure, with which you can:
% - draw around a set of sensors
% - click in the drawn box to produce a plot of the time series
% - on the time series plot you can draw a time window
% - and click in the window to create a topoplot averaged over that time
% window (which is itself interactive....!)

S2=[];
S2.oat=oat;
S2.stats_fname=oat.first_level.results_fnames{1};
S2.modality='MEGPLANAR';
S2.first_level_contrast=[3];
S2.cfg.colorbar='yes';
S2.cfg.zlim = [-5 5];
S2.view_cope=0;

% calculate t-stat using contrast of absolute value of parameter estimates
[cfg, data]=osl_stats_multiplotTFR(S2);
title([oat.first_level.contrast_name{S2.first_level_contrast}]);


%% CREATE A TOPOPLOT AVERAGED WITHIN A TF WINDOW
%
% This section calls another fieldtrip function which creates a topoplot 
% averaged over 130 to 160 ms, and 8 to 12 Hz.

cfg.xlim        = [0.13 0.16]; % time window in secs
cfg.ylim        = [8 12]; % freq window in Hz
cfg.interactive = 'no';
figure; ft_topoplotTFR(cfg,data);
title([oat.first_level.contrast_name{S2.first_level_contrast}]);


%% CHALLENGE!
%
% Can you modify a copy of this script to run an OAT analysis to look for
% changes within a single frequency band? Specifically, set up an analysis
% which looks for changes in power between 5 and 20Hz over the course of
% the epoch.
%
% The answer can be found in osl_example_sensorspace_oat_tf_answer.m
